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ABSTRACT 


Element-by-element approximate factorization, implicit-explicit and adaptive 
implicit-explicit approximation procedures are presented for the finite-element 
formulations of large-scale fluid dynamics problems. The element-by-element 
approximation scheme totally eliminates the need for formation, storage and inversion 
of large global matrices. Implicit-explicit schemes, which are approximations to 
implicit schemes, substantially reduce the computational burden associated with large 
global matrices. In the adaptive implicit-explicit scheme, the implicit elements are 
selected dynamically based on element level stability and accuracy considerations. 
This scheme provides implicit refinement where it is needed. 

The methods are applied to various problems governed by the 
convection-diffusion and incompressible Navier-Stokes equations. In all cases 
studied, the results obtained are indistinguishable from those obtained by the implicit 
formulations. 
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CHAPTER 1 


INTRODUCTION 

Significant improvements have been made in computer memory and speed in 
recent years. However, the existing computer power is still far behind the demands 
from scientists for large-scale fluid dynamics calculations. Implicit schemes, which 
are desirable for their stability and accuracy properties, lead to large global coefficient 
matrices which need to be formed, stored and inverted. Direct management of such 
matrices is very difficult for large-scale, geometrically complicated two-dimensional 
problems and is virtually impossible for large-scale three-dimensional problems. For 
example, a cubic domain with 30 nodal points along each side and with four 
unknowns at every node will produce some 108,000 equations with 7.8 x 10 8 
entries in the global coefficient matrix. Even with the currently available 
supercomputers it is almost hopeless to try to handle directly such a large equation 
system. 

To overcome the shortage of computer power, many researchers have developed 
algorithms for large-scale problems. In their application of domain decomposition 
methods, Glowinski, Dinh and Periaux [1] successfully coupled the incompressible 
viscous flow and incompressible potential flow models employed in different 
subdomains. A conjugate-gradient method, which is basically suitable for the 
symmetric and positive-definite systems, was employed to solve the variational 
problem. For the nonsymmetric and nonpositive-definite systems, one has to find an 
appropriate preconditioner for each problem; this needs much sophisticated work 
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[2-4]. A three-dimensional flow simulation with 128 3 nodal points was made by 
Rogallo [5] with the alternating-direction method. A new class of algorithms in 
numerical linear algebra which take advantage of the parallel computation capabilities 
of modem computers also provides hope [63- 

In this report, we present element-by-element (EBE) approximate factorization, 
implicit-explicit (IMEX), and adaptive implicit-explicit (AIE) schemes for large-scale 
computations in fluid dynamics. Compared to implicit methods, these schemes 
possess, to a great extent, similar desirable stability and accuracy properties, yet result 
in substantial reduction in computer memory and CPU time demands. 

The element-by-element approximate factorization is in some sense related to 
domain-decomposition [1,7] and alternating-direction ( or operator-splitting ) [8-11] 
schemes. These schemes, contrary to the EBE method, are based on global 
approximations. The EBE method was first proposed by Hughes, Levit and Winget 
[12,13] with applications to transient heat conduction, structural and solid mechanics 
problems. Preliminary application to fluid mechanics problems within the context of 
the compressible Euler equations were presented in Hughes, Levit, Winget and 
Tezduyar [14], In this report, we present our EBE formulation for problems with 
nonsymmetric, nonpositive-definite spatial differential operators. Currently we focus 
on the convection-diffusion and the incompressible Navier-Stokes equations. In the 
EBE formulation, global coefficient matrices are approximated by sequential product 
of much simpler matrices which are based on the most natural unit in the finite-element 
method, the element-level matrix. Every calculation is done at the element level. The 
method keeps the versatility of the finite-element formulation in its easy adjustment to 
irregular meshes. Other advantages of the finite-element method such as easy 
implementation of the boundary conditions and the source terms are also retained in the 
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EBE formulation. 

Our implicit-explicit approximation schemes for fluid dynamics problems are 
based on the work of Hughes and Liu [15,16] which was for solid mechanics and heat 
transfer applications. We consider problems governed by the convection-diffusion 
and the incompressible Navier-Stokes equations. In this approach, for the elements 
which are designated to be implicit, the element level matrices are kept as they are, 
whereas for the explicit elements the element level matrices are approximated by a 
diagonal matrix. The implicit element-explicit element decision is based on the local 
stability criterion. 

In the adaptive implicit-explicit scheme the implicit element-explicit element 
decision is made dynamically. The stability and accuracy criteria applied to each 
element, based on the solution from the previous iteration or previous time step, 
determine whether an element needs to be implicit In this approach we can adaptively 
have implicit refinement where it is needed for stability and accuracy. Compared to 
other adaptive schemes which are based on element redistribution or element 
subdividing [17], the AIE approach involves minimal bookkeeping and no geometric 
constraints; therefore the method is very easy to implement 

In Chapter 2, model problems are described and their spatial and temporal 
discretizations are given. In Chapter 3, the stability and accuracy analysis needed for 
the development of the IMEX and AIE schemes is performed. The EBE, IMEX, and 
AIE schemes are discussed in Chapters 4 and 5. Numerical results are presented in 
Chapter 6 and the conclusions are given in Chapter 7. 
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CHAPTER 2 


MODEL PROBLEMS - SPATIAL AND TEMPORAL 
DISCRETIZATIONS 

We consider model problems governed by the convection-diffusion equation 
and the vorticity stream-function formulation of the two-dimensional incompressible 
Navier-Stokes equations. A formal statement of the problems and the corresponding 
spatial and temporal discretization for both cases are given below. 

Convection-Diffusion Equation 

Let ft and T denote the computational domain and its boundary. 
Time-dependent convection-diffusion of an unknown scalar function d> is governed by 
the following set of differential equations, boundary conditions, and initial condition: 

<D, t + u-V $> = V-(KV<I>) + f, on ftx]0,T[ (2.1) 

<D(x,t) = g(x,t), Vxe r g ,teJ0,T[ (2.2) 

n • k VO (x,t) = h(x,t), V x € T h , t € ] 0,T [ (2.3) 

and 

O ( x , 0 ) « O 0 ( x ), on ft (2.4) 

where the velocity field u » u ( x ) is given and k is the conductivity matrix. The 
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source term is given asf = f(x,t). The outward unit vector normal to the boundary 
r is denoted by n, whereas g, h and <I> 0 are prescribed functions. r g and are the 

mutually exclusive but complementary subsets of the boundary T with Dirichelet and 
Neumann boundary conditions respectively. That is, 

0*r g nT h (2.5) 

r=r g ur h ( 2 . 6 ) 

Multiplying equation (2.1) by the weighting function w and applying Green's 
theorem, we obtain the following weak form. 


w<X>, t dQ+ wu*V<X>dn + Vw 

Jn Jo Jn 

= w h dT + wfdfl 

Jr h j£i 


( K VO) d£2 


(2.7) 


and 


h 

Ja 


w ( O - Oq) dfi s 0 


( 2 . 8 ) 


Finite-element spatial discretization of equations (2.7) and (2.8) leads to the following 
semi-discrete equations: 


M0>+ C®= F 

(2.9) 

O(0) = 4> 0 , 

(2.10) 
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where M, C and F are the "mass" matrix, "stiffness" matrix, and the generalized 

• 

"force" vector , respectively; O and <X> represent the dependent variable and 
its temporal derivative at the nodes. The information regarding the initial condition is 
contained in vector <X> Q . 

Burger's equation, given below, can be treated as a special case of the 
convection-diffusion equation. 


<& n +OO, x = 0 V x 6 ] 0,L [, te ] 0,T[ (2.11) 

This simple nonlinear hyperbolic equation is used to study the numerical performance 
related to shock formulation and entropy condition. A finite-element discretization of 
this equation leads to the following semi-discrete nonlinear equation system, 

M0 + N(0) = F, (2.12) 

where N ( <1>) is a nonlinear function of O. 

Vorticitv Stream-Function Form of the Two-Dimensional 

Incompressible Navicr-Stokcs Equations 

The incompressibility condition and the nonlinear convection term constitute 
some of the major difficulties for the numerical methods associated with the 
Navier-Stokes equations. By employing a vorticity stream-function formulation, one 
can easily handle the incompressibility condition. However, the extension of the 
numerical algorithm to three-dimensional problems is rather cumbersome and the 
treatment of the boundary conditions on no-slip surfaces and at out-flow boundaries is 
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complicated ( see, e.g., [18] ). The field equations consist of a time-dependent 
transport equation for the unknown vorticity function 0) and an equation which relates 
the unknown stream function to vorticity. They are given as follows: 


(0 n + u-Vco = v V 2 (d on £2 x ] 0,T [ (2.13) 


and 


V 2 = - 0), on Q (2.14) 

where v is the kinematic viscosity. The velocity field u, which is now an unknown, is 
related to the stream function by the following equations: 

Uj = 8'F/8x 2 (2.15) 


and 

u^= - d*¥ /8xj 


(2.16) 


Note that equations (2.15) and (2.16) lead to a zero-divergence condition. 


V-u = 8 2 'P /d *2 - i) 2 '? /dxj dx 2 = 0 


(2.17) 


In other words, the incompressibility condition on u is assured by the definition of 'P. 
The boundary conditions for (0 and are given below: 


co(x,t)= g(x,t). 

V X € 

r ~ 
1 g’ 

te]0,T[ 

(2.18) 

v n • Vo) ( x,t ) = h ( x,t ), 

V X € 


t e ] 0,T [ 

(2.19) 

'V ( x,t ) = g ( x,t ), 

V X € 

r g* 

te ]0,T[ 

(2.20) 
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n-V’F ( x,t ) = h ( x,t ), 


Vxer h , t € ] 0,T [ (2.21) 


where and Pf are the mutually exclusive but complementary subsets of the 

boundary T with Dirichelet and Neumann boundary conditions for co, whereas T and 

© 

Tjj are similar boundary subsets for that is 


0 = r~nr£ 

(2.22) 

r-rjurj; 

and 

(2.23) 

0 " r g nr h 

(2.24) 

r - r g^ r h 

(2.25) 


The difficulty associated with the lack of boundary conditions for the vorticity and its 
normal derivative on no-slip surfaces can be handled by a discrete Green's formula 
approach which is described in Tezduyar, Glowinski, and Glaisner [19]. 

The initial condition for the vorticity is given as 

o>(x,0) = cd 0 on ft (2.26) 

Multiplying (2.13) and (2.14) by weighting functions w and w respectively, the 
following weak forms can be obtained via the Green's theorem: 



r 

w h dT 

Jnr 


~ I ~ 

w u-Vco dQ + Vw • v V to di2 

a Jn 


(2.27) 
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and 


wcodfi + 

VcoV'FdQ = 

r 

Jn 

Ja 

*/ • 


w h dT 


(2.28) 


Finite-element spatial discretization of equations (2.27) and (2.28) leads to the 
following set of equations: 


M to + K to = F, 


(2.29) 


and 


to ( 0 ) = to 0 . 


-Mco + K'F = F, 


(2.30) 

(2.31) 


where M, K, F, M, K, and F are the matrices and vectors resulting from the spatial 
discretization; co, to, and 'F represent the vorticity, its temporal derivative and the 
stream function at the nodes. The initial condition for to is represented by co 0 . 


Spatial Discretization 

Due to the existence of nonsymmetric spatial operators ( convection terms) in 
our model problems, the "best approximation" property of the Galerkin finite element 
formulations [20] is lost, and therefore spurious oscillations can be encountered in 
convection-dominated problems. Streamline-Upwind/Petrov-Galerkin (SUPG) 
formulations are employed for the spatial discretization of our model problems. These 
formulations are well-known for their robust and accurate performance for problems 
with nonsymmetric spatial operators. Several papers can be found in scientific 
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literature on this type of schemes ( see, e.g., [21, 22]). For the purpose of this 
thesis, we briefly describe the SUPG scheme employed in terms of the weighting 


function basis N a ; this description is given by 
0*0 

N, = N a + °2t h /2 *•' VN t (2.32) 

or 

N m = N a + x u*VN a , (2.33) 

where the subscript "a" refers to an element node, N a is the solution function basis, h 
is the "element length" in the direction of u , and s is the unit vector in u direction. 

The "algorithmic Courant number" C 2x and the "algorithmic time constant" x are 
related by the following expression, 

C^-Hu Il2x/h. (2.34) 



where the vector a is the temporal derivative of the vector v. N (v) is a vector- valued 
function of v which we will, in general, assume to be nonlinear. 

The tangent " stiffness matrix" K is defined as 
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K = d N (v) / 9 v (2.37) 

A predictor/multi-corrector transient integration algorithm [23] is employed to 
solve equations (2.35) and (2.36). Let a subscript denote the time step and a 
superscript denote the iteration step. The algorithm can be summarized as follows: 


1 . given v 0 find a 0 from 

M a 0 + N (Vq) * F 0 

n — » n+1 

2. predictor stage 

v n + i 0 ssV n + O-«)^ta n 

Vi° = 0 

i-» i+1 

3. corrector stage 

M* Aa„,i = R n .i' 

Vi l+1 = Vi i+ 4 Vi' 

and 

Vi i+1 = W + a AtAa^j*. 


(2.38) 


(2.39) 

(2.40) 


(2.41) 

(2.42) 

(2.43) 

(2.44) 


where a e [ 0,1] is a parameter which controls the stability and accuracy. 
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Convergence is checked by inspecting the norm of R n+1 i and Aa n+1 ‘ . 


Remark : 

1 . A consistent derivation for M* gives the following expression, 

M* = M + a At K (2.45) 

If left as it is, this expression for M* leads to an implicit formulation which 
requires only one correction for linear systems. However, for nonlinear 
systems, one needs to have as many corrections as the convergence criterion 
dictates. 

2. A choice of 

M*=M U (2.46) 

where M L is a lumped version of matrix M , leads to an ex plicit formulation. 

Explicit formulations in general are less stable, less accurate, but also less costly 
(less memory and less CPU time) compared to implicit schemes. The 

conditional stability of these methods are usually expressed in terms of a limit 

on the element Courant number. The element Courant number C At is defined 
as 

C^= llu IIAt/h (2.47) 

3. In the SUPG formulation, the choice of 

T = aAt (2.48) 

leads to a symmetric positive-definite M* for pure convection problems. 
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Proof : 

A general convection equation is given as 

<D, t + u-V <D= f (2.49) 

The following semi-discrete ( spatially continuous) formulation can be obtained 
from equation (2.49), 

( °n + l • 1 At + uV < (J-tt) ^n+i) 

= (l-a)f n + af n+1 (2.50) 

For the weak form of equation (2.50), we employ a weighting function 
w which is given as 

w =w + xu*Vw, (2.51) 

where w and <J> come from the same space. Note that for this pure convection 
system, a Dirichlet type boundary condition for <I> is required on the part of 
the boundary where the information comes from outside the domain, and w has 
to satisfy the homogeneous form of the same boundary condition. That is, 

O ( x ) = g ( x ) Vxe {xlxeT, n(x)-u(x)<0} (2.52) 

and 

w ( x) = 0 Vxe (xixe r,n(x)-u(x)<0 ) (2.53) 

The weak form of equation (2.50) with the weighting function of equation 
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(2.51) can be written as 


(w, ^n+i ) + a At ( w > u V ) + a At (u-Vw, O n+1 ) + 

( a At ) 2 (u-Vw, u-V <|> n+1 ) = right hand side (2.54) 

where the bilinear form ( • , • ) is defined as 


( w,<D) 


f w <X> dQ 


Jn 

Note that only the left hand side of equation (2.54) corresponds to M*. 
Rewriting equation (2.54), we get 

(w + a At u-Vw, <D n+1 + a At u-V<X> n+1 ) = righe hand side 


(2.55) 


(2.56) 


We need to show that the left hand side of equation (2.56) has a symmetric 
positive-definite form. It is obvious that this bilinear form is symmetric and that 


(w + aAtu-Vw, w + aAtu*Vw) £ 0 

(2.57) 

We have to show that if 


(1/ a At ) w + u-Vw * 0, 

(2.58) 

then 


w = 0 

(2.59) 
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Note that equation (2.58) is nothing more than a linear, steady-state, 
convection-reactiontypc equation for w. Since (1/aAt ) is always positive, 
the reaction is a consumption on w. If the value of w at any point can be traced 
back to a boundary point, because of equation (2.53), w = 0; if not, w = 0 
because of consumption. This completes our proof, based partly on physical 
reasoning, that M* is symmetric positive-definite. 
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CHAPTER 3 


STABILITY AND ACCURACY ANALYSIS 

In convection-dominated problems, Galerldn formulation, just like the classical 
central-difference schemes, generates numerical oscillations especially in the presence 
of sharp layers. To minimize these oscillations, one can refine the grid at the 
oscillation zones, but then must pay for increased memory and CPU costs. Another 
option is the classical upwind scheme which leads to forward- or backward-difference 
treatment of the convection terms depending on the direction of the convection 
velocity. With this option, while the wiggles are minimized, usually an excessive 
amount of artificial diffusion is introduced into the numerical simulation. To explain 
this, let us consider a simple one-dimensional convection problem, 

<J>, t +ud>, x = 0, (3.1) 

where u > 0 is the convection velocity. For the node numbered i at time step n+1, the 
upwind scheme leads to the following fully discrete equation, 

(1/At) (<D i n+1 - <& ") + u (1/h) (<&."♦! - O. j»+i) = 0, (3.2) 

in which At is the time integration step and h is the spatial discretization step. It is 
clear that equation (3.2) can be rewritten as 
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(1/At) (O.n+1 - 0.») + u (l/2h) (O i+1 ”+i - 0. ^+1) 


- (u h / 2) (1/h 2 ) ( <D i+1 n+1 + 0. ^+1 - 2<J>. n+1 ) = 0 (3.3) 

Returning to the continuous counterpart of equation (3.3), we obtain 

O n + uO^- (uh/2)O, xx = 0 (3.4) 

It can be observed that, compared to equation (3.1), an artificial diffusion term with 
coefficient value of (uh / 2) has been introduced into the continuous equation by the 
upwind scheme. For high convection velocities and large elements, the artificial 
diffusion becomes significant This results in simulations which depart substantially 
from the real physical phenomenon. 

Streamline-Upwind/Petrov-Galerkin formulations [21,22] keep the numerical 
oscillations and the level of artificial diffusion to a minimum. Weighting functions 
leading to the SUPG formulations are described by equations (2.1)-(2.3). Tezduyar 
and Ganjoo [21] have investigated the various choices of the term in the SUPG 

formulations by a phase accuracy and amplitude analysis for a one-dimensional 
transient pure convection equation (3.1). These choices of the term are 


C 2t= 1 

(choice 1) 

(3.5) 


(choice 2) 

(3.6) 

0^= 2/Vl5 + (l- 2/Vl5)C At 

(choice 3) 

(3.7) 


17 


In [21], the exact solution of equation (3.1) is assumed to be of the form, 

0(x,t) = e^ +i ®> t e- aa (3.8) 

and the finite-element solution is of the form, 

<D(xi,tn) = e' 3 ^, (3.9) 

where q and q are the exact and approximate damping rates; co and go represent the 
exact and approximate frequencies. The wave number is denoted byk(k = 27t/X, X 
is the wave length ). The significant quantities in this analysis are the algorithmic 
damping ratio (ADR = ^ / co ) and the frequency ratio (FR = S / go ). A good scheme 
should have an ADR approaching zero and a FR approaching unity. From the rates 
these quantities approach these limits, one can infer the order of accuracy of the 
algorithm. 

Let h denote the element length; a dimensionless wave number is defined as 
q = kh (3.10) 

This is a measure of the spatial refinement of the numerical method. Figures 3. 1-3.3, 
adopted from [21], show the ADR and FR values versus q for a set of element 

Courant numbers and various choices of C^. The order of accuracy can be 

revealed from the ADR and FR curves. A scheme is first-order accurate if either of the 
curves has a finite slope as q approaches zero. It is at least second-order accurate if 
both curves have zero slopes as q approaches zero. It can be seen that all implicit 
schemes have at least second-order accuracy. In fact, choice 3 leads to a fourth-order 
accurate scheme [24] as C* approaches zero. For explicit one-pass schemes, 
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SUPG IMPLICIT 

c o - 2/vi 6 - (i-2/v15)C 

♦ 



DlMEN Si ON LESS WAVE NUMBER (q) 


SUPG IMPLICIT 

f C ^ - 2/Vi 5 + (1-2/vi5)C a 



Fig. 3.1 Stability and accuracy characteristics for the implicit SUPG formulation. 
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Fig. 3.2 Stability and accuracy characteristics for the explicit one-pass SUPG 
formulation. 
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choice 2 gives second-order accuracy. All explicit two-pass schemes are second-order 
accurate. 

Based on a single degree of freedom model analysis performed by Tezduyar 
and Hughes [25] for the one-dimensional convection equation, we obtain the 
following stability conditions: 


1. Implicit scheme: 

All choices of are unconditionally stable for a £ 1/2 and r £ 1/4 where r is a 

parameter which depends on the numerical integration technique used for the time 
dependent terms. Typical values are r = 0 for lumped mass, r = 1/6 for consistent 
mass and r = 1/4 for PG (pade) mass ( see [21] ). 


2. Explicit one-pass schemes: 

All choices of are conditionally stable; the stability criterion yields 


^'At < 


2C* 


2-vd-C*) 


(3.11) 


where v = 1- cos q. Note that, for Galerkin schemes ( = 0 ), the method is 

unconditionally unstable. It can be shown that for all choices of C 2t the stability 
limit is 


^<1 

3. Explicit two-pass schemes: 

The following inequality must be satisfied for stability: 


(3.12) 
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2 ( A r Bj - Aj B 2 ) + (A^ + Aj2 ) ( B^ + B 2 2) < 0, (3.13) 

where 

Ai = C At q T v (3.14) 

AjsC^w (3.15) 

B l = aC At C 2t v * 1_2rv (3.16) 

and 

B 2 » cx C At w - w/2, (3. 17) 


in which w = sin q. The worst case occurs for q = 7 t (which leads to v = 2 and w = 
0); then equation (3.13) implies 

4C*C 2t (2aC At C 2t -l-4r)x 

{ 2 a [ C At - (1+ 4r)/4a] 2 + 1 - (l+4r) 2 /8a } < 0 (3.18) 

It is obvious that under the choices of a £ 1/2 and r £ 1/4, 

1 -(l+4r) 2 /8a £0, (3.19) 


then equation (3.18) leads to the following stability condition, 

1 +4r 

— (3.20) 

2a 

For various choices of C^, the stability limits can be expressed as follows: 
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choice 1: C At < — * + ■ (3.21) 

2a 

choice 2: C At < ( — - + - - ) 1/2 (3.22) 

2a 

-2/Vl5 + V4/15 + 4 ( 1- 2/V15) ( l+4r/2a) 

choice 3: C At < (3.23) 

2(l-2/Vl5) 

Similar analysis has been done in one-dimensional diffusion equation for the Galerkin 
finite-element formulation. The conditional stability of the explicit schemes are 
expressed in terms of a limit on the diffusion Courant number. The diffusion Courant 

number C K is defined as 

C K = 2 k At / h 1 2 , (3.24) 

where k is the diffusion coefficient. The stability and accuracy properties can be 
summarized as follows: 

1. Implicit scheme: 

a) Unconditional stability is assured for a £ 1/2 and r £ 1/4. 

b) The algorithm is second-order accurate. It becomes fourth-order accurate with 

the choice of r * 1/12. 
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2. Explicit one-pass scheme: 

a) It is conditionally stable. The stability criterion is 

Ck* 1 

b) The algorithm is first-order accurate. 


3. Explicit two-pass scheme: 

a) It is conditionally stable . The stability criterion is 




1+4 r 
2a 


b) The accuracy is of the same order as the implicit scheme. 
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(3.25) 


(3.26) 


CHAPTER 4 


ELEMENT-BY-ELEMENT (EBE) APPROXIMATE 
FACTORIZATION 

The equation system of (2.42) can be rewritten as follows, 

A x = b, (4.1) 

where x is the increment vector Aa, b is the residual vector R, and A is the coefficient 
matrix M*. 

A parabolic regularization of equation (4. 1) can be expressed as 

W dy / d0 + A y = b, (4.2) 

where 0 is a dimensionless "pseudo-time" and 

y ( 0 ) = 0 (4.3) 

We assume that 

x = lim y ( 0 ) (4.4) 

e_»oo 

The choice of W depends on the properties of A. For a symmetric 
positive-definite A, Hughes, Levit and Winget [12,13] proposed W to be the 
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diagonal part of A, 

W = diagA (4.5) 

This was found to be effective for all problems studied in [12,13]. Alternately W can 
be chosen to be the lumped mass matrix mentioned in Chapter 2. That is, 

W = M l (4.6) 

This choice was proposed by Hughes, Winget, Levit and Tezduyar [14]. 

We employ an unconditionally stable pseudo-time stepping algorithm given 
below: 


(W + aA8A)Ay m - A9r m 

(4.7) 

£ 

< 

1 

x> 

II 

£ 

(4.8) 

y m+ i = y m +A ym’ 

(4.9) 


in which a e [1/2, 1], The scheme is backward-Euler with a = 1, and 
Crank-Nicolson with a = 0.5. Equation (4.7) can also be written as 

Ay m = D ( I + a DAD ) > D r„, (4.10) 

when I is the identity matrix and 

D = ( W* 1 A0 ) iri (4.11) 

Note that equation (4.11) requires W to be positive-definite. While the choice of 
equation (4.5) does not guarantee this when M* is not positive-definite, the alternate 
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choice of equation (4.6) does. One also needs to note that for problems in fluid 
mechanics, M* in general is not symmetric, positive-definite; however, as mentioned 
in Chapter 2, the choice of equation (2.50) assures that it is. 

Since we are interested in the steady-state solution of equation (4.2), usually the 
backward-Euler scheme is adopted to avoid undamped and oscillating solutions as the 
condition number of the problem deteriorates or as large pseudo-time steps are chosen 
[26], 

Rewriting equation (4. 10) with the choice of a * 1, we obtain 

Ay m = D (1+ DAD ) 4 D r m (4.12) 

The EBE approximate factorization is based on the following expressions: 

1. one-pass EBE approximation: 

net 

( I + DAD) - n ( I + DA e D ), (4.13) 

e*l 

in which nel is the total number of elements in the domain, and A e is the 
contribution of e* element ( or subdomain) to A. Clearly, this approximation 
depends on the element ordering. 

2. two-pass EBE approximation: 

nel 1 

(I + DAD) - 11(1+1/2 DA e D) II (1+ 1/2 DA e D) (4.14) 

e=l e=nel 
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It should be noted that the approximation of the one-pass EBG deteriorates as a 

large A6 is taken to assure a fast steady-state solution of equation (4.2). The 
two-pass EBE approximation has a better performance in this situation. 

We update y m+1 according to the following formula, 

y m+l = y m + S A y m’ ( 4 - 15 ) 

where s is a search parameter obtained by minimizing II r m+1 II 2 with respect to s, 

s =(AAy m )T m /IIAAy m ll 2 (4.16) 


R e ma rk: 

1. The need for the formation and storage of the global matrix A has been eliminated. 
There is no need to store the element level matrices A e 's; however, if desired, 
instead of recomputing for each EBE iteration, the element level matrices can be 
stored. Even then, the storage requirement is far less than that of a global matrix. 
For instance, in a cubic grid with N nodes at each side and one unknown at each 
node, 2N 5 entries must be stored in the global matrix for the implicit calculation. 
EBE needs only 64 entries in each element matrix. Even if all element level 
matrices are stored, the total number of entries is only 64N 3 . 

2. The EBE approximate factorization procedure is parallelizable. This aspect 
makes EBE favorable especially since parallel computations are expected to 
play a more and more important role in computational fluid dynamics. 

3. It must be well understood that if the EBE procedure converges, it converges 
to the solution of equation (4.1), which is the equation system of the implicit 
method. 
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CHAPTERS 


IMPLICIT-EXPLICIT APPROXIMATION SCHEMES 

In solving partial differential equations via finite-element method, since the 
solution may vary significantly in the domain, a nonuniform mesh is usually employed 
for stable, accurate, and efficient solution procedures. A finite element analyst can 
vary the element sizes in his/her nonuniform mesh anyway he / she pleases but usually 
according to some criteria based on past experience, physical intuition, and 
mathematical insight Implicit schemes, though stable and accurate, may cost too much 
(CPU time and storage) to compute. To guarantee stable solutions for explicit 
schemes, one has to limit the time discretization step size based on the smallest element 
size in the entire domain. This may place a rather prohibitive constraint on the 
explicit integration scheme. The implicit-explicit approximation schemes provide an 
effective way to overcome these drawbacks. The basic idea associated with the 
implicit-explicit schemes can be well demonstrated by the following simple 
one-dimensional initial /boundary-value problem: 

O h +u^-k^ x =0, V xe ] 0,1 [, t e [ 0,T ] (5.1) 

0>(0,t)=l, V 1 6 [ 0,T ] (5.2) 

O(l,t)»0, V t e [ 0,T ] (5.3) 

and 

0(x,0) = 0. Vxe]0,l[ (5.4) 
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The steady-state analytical solution to equation (5.1) is 


<D = 


gUX/K. gU/K 


1 -e"^ 


(5.5) 


For sufficiently large u/k, the solution is essentially unity in the interior of the domain, 
decaying through a boundary layer near x = 1 to assume the boundary value ofd> = 0. 
In our computation, a finer mesh is employed in the right region of the domain as 
shown in Figure 5.1(a). A constant time step size is chosen such that the element 
Courant number is 0.7 in the left and 1.4 in the right regions. Figure 5.1(b) shows the 
solution obtained by the implicit scheme. Figure 5.1(c) shows the unstable solution 
obtained by the explicit one-pass scheme; this is due to the high Courant number in the 
right region. In the implicit-explicit scheme, we treat the elements in the right region 
implicitly and the ones in the left region explicitly. The solution is stable and accurate 
as shown in Figure 5.1(d). 

Let £ be the set of all elements, e = l,2,...,nel. The assembly of the global 
matrix M* can be expressed as follows. 


M*=A (m*) e , (5.6) 

eeC 

where A is the assembly operator. 

Let £j and £ E be the subsets of £ corresponding to "implicit elements" and 
"explicit elements" respectively, such that 


£ - £ju£g 

and 

0 = £j n £ e 


(5.7) 

(5.8) 
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Fig. 5.1 One-dimensional boundary layer problem: (a) Finite element mesh. 
Solutions by various schemes: (b) Implicit scheme, (c) Explicit one-pass scheme, 
(d) Implicit-explicit one-pass scheme. 
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Rewriting equation (5.6), we get 


M* = A (m*) e + A (m*) e (5.9) 

ee 8j ee Eg 

Implicit-explicit approximation [15,16] is based on the replacement of (m*) e by the 
element level lumped mass matrix (m) e L , 

M*= A (m*) e + A (m) e L (5.10) 

cs 8j ee Eg 

Let M 1 and K 1 denote respectively the mass matrix and stiffness matrix obtained 
by assembling the implicit elements and M E denote the diagonal mass matrix from the 

explicit elements ( aAt K E is neglected). An equivalent form of equation (5.10) can be 
expressed as 

M* = M e + M 1 + aAt K 1 (5.1 1) 

Many finite element simulations involve movement of fronts in the computational 
domain. The unknown variable can be the temperature in the heat equation, the 
vorticity in the Navier-Stokes equations, or the oil concentration in the enhanced oil 
recovery processes. Due to the presence of the convection operator and sharp layers in 
the solution, the solution is prone to instability and inaccuracy around the moving 
fronts. In such cases implicit-explicit schemes can be employed by concentrating the 
implicit elements around the moving fronts. The implicit elements need to somehow 
follow the moving fronts as shown in Figure 5.2. 
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Fig. 5.2 Implicit element distribution in the moving front problem. 
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We propose that the sets £j and £g are determined dynamically. The 
determination will be based on several criteria including stability and accuracy 
considerations for all types of transport phenomena present Stability and accuracy 
characteristics of algorithms are described not only in terms of the Courant numbers 

(C At and C K ) but also in terms of the dimensionless wave number ( see Figures 

l(a)-3(f) ). For the stability consideration, the elements for which the Courant 
numbers exceed the stability limit of the explicit schemes as described in Chapter 3 arc 
assigned to be implicit while explicit elements are introduced elsewhere to save 
computer memory. The accuracy criterion is determined by the dimensionless wave 
number of the solution from the previous time step (or iteration). We can achieve this 
by defining a determination criterion which is based on some measure of the jump in 
the solution or jump in the flux across an element under consideration. The jump 
values are computed based on the previous time step (or iteration). Currently we 
employ the following definition for the jump in the solution across an element, 

[ O ] = max ( <D a ) - min ( <I> a ), (5.12) 

a a 

where "a" is the element node number. Note that the threshold value for this jump 
which makes the element eligible for the implicit set £j should be based on the global 
scaling of the solution field. Such a global scaling idea is closely related to the global 
scaling concept of the "discontinuity capturing" scheme described in Tezduyar and 
Park [27]. 

In the adaptive implicit-explicit (AIE) approach, one can have a high degree of 
refinement throughout the mesh but can raise the implicit flag only for those elements 
which are proposed to be treated implicidy. 
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Compared to other adaptive concepts such as adaptive element redistribution or 
adaptive element-subdividing, the AIE scheme is far easier to implement because it 
involves no geometric changes; the bookkeeping involved is minimal. 

Remark : 

1. By appropriately numbering the nodes based on the distribution of the implicit 
elements, one can substantially reduce the bandwidth of the global matrices. 

2. By the aid of bandwidth optimizer packages available in the market, one can 
renumber the nodes every time step ( or every iteration) in the AIE calculations 
to efficiently save the computer memory. 
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CHAPTER 6 


NUMERICAL EXAMPLES 

The EBE, implicit-explicit, and AIE schemes have been tested on various 
problems governed by the model equations stated in Chapter 2. The results are 
compared with those obtained by the implicit and explicit schemes. The implicit 
solutions were shown by Ganjoo [28] and Glaisner [29] to be in agreement with the 
results from various past publications. 

One-Dimensional Advection of a Cosine Wave 

This problem consists of pure advection of a cosine wave from left to right A 
cosine wave of unit amplitude is initially set in the left part of the domain. A uniform 
mesh containing SO elements and 102 nodal points is employed. A homogeneous 
Dirichlet boundary condition is specified on the extreme left, and a homogeneous 
Neumann boundary condition is imposed on the extreme right The advection velocity 
is 1.0 and the time step is such that the element Courant number is 0.6. Figures 
6.1(a)-6.1(d) compare at the end of the same time interval the results obtained by 
various schemes. Figure 6.2 shows the distribution of implicit elements for the AIE 
calculation at various time steps. The amplitude of the cosine wave decays to 0.969 
for the implicit and EBE schemes, and to 0.965 for the adaptive implicit-explicit 
one-pass (AIE-1) scheme. A considerable decay of amplitude (to 0.769) is observed 
for the explicit one-pass (EXP-1) scheme. 
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Fig. 6.1 One-dimensional advection of a cosine wave: Solution obtained by various 
schemes at t= 0.0 and 0.408 (a) Implicit scheme, (b) EXP-1, (c) AIE-1. (d) EBE. 
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Fig. 6.2 One-dimensional advection of a cosine wave: Implicit elements in the ATF. -1 
calculation at various time steps. 
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A discontinuity is advectcd from left to right under the same set-up conditions as 
the previous case. Figures 6.3(a)-6.3(d) show the results from various schemes at the 
end of the same time interval. Similar overshoots are found in the implicit, EBE, and 
AIE-1 schemes. Undershoot is observed in the implicit and EBE schemes. By 
changing the jump criterion described in Chapter 5, the same undershoot can be 
obtained in the AIE-1 scheme. The discontinuity becomes smeared in the EXP-1 
scheme. Figure 6.4 shows the distribution of implicit elements for the AIE calculation 
at various time steps. 

Two-Dimensional Advection of a Cosine Hill 
(Translating Puff) 

Two-dimensional advection and rigid body rotation of a cosine hill are typical 
examples for testing algorithms in convection-dominated problems ( see, e.g„ [30] ). 
The translating puff problem consists of advection of a cosine hill from the extreme 
left to the right. Two meshes are tested: a uniform mesh with 30 x 30 elements in a 
lxl domain and a nonuniform mesh with 45 x 30 elements in a 1 x 0.75 domain. For 
the uniform mesh an initial cosine hill profile with unit peak amplitude and base radius 

of 0.2 is centered at (x^Xj) = (0.267, 0.5). The diffusion coefficient is set to be 10"*; 

the advection velocity is unity in Xj-direction, and the time step is adjusted to give a 
Courant number of 0.6. An homogeneous Dirichlet boundary condition is specified 
on all boundaries except at Xj = 1.0 where an homogeneous Neumann boundary 

condition is imposed. Figures 6.5(a)-6.5(c) show the results at various times obtained 
by AIE-1. Figures 6.6(a) - 6.6(c) show the distribution of the implicit elements at the 
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Fig. 6.3 One-dimensional advection of a discontinuity: Solution obtained by various 
schemes at t= 0.0 and 0.408 (a) Implicit scheme, (b) EXP-1, (c) AIE-1. (d ) E BE. 
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Fig. 6.4 One-dimensional advection of a discontinuity: Implicit elements in the ATF.-1 
calculation at various time steps. 
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Fig. 6.5 Elevation plots for the translating 
puff on a uniform mesh obtained by the 
AIE-1 scheme.(a) Initial condition. 

(b) at t = 0.3. (c) at t = 0.72 . 


Fig. 6.6 Distribution of the implicit 
elements for the A EE- 1 calculations 
of the translating puff on a uniform 
mesh, (a) at t = 0. (b) at t = 0.3. 

(c) at t = 0.72. 
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corresponding times for the AIE-1 computations. Figures 6.7(a)-6.7(b) show the 
results at various time steps obtained by EXP-1. The peak amplitude value of the 
cosine hill at t = 0.72 is 0.972 for the implicit and EBE schemes and 0.974 for AIE-1. 
EXP-1 gives a stable solution with poor accuracy. 

For the nonuniform mesh, the element length in the left region of the domain is 
half of that in the right region. The time step is chosen such that the Courant number 

is 1.8 in the left and 0.9 in the right The initial profile is centered at (Xj.Xj) = 
( 0.233, 0.375 ) with base radius of 0.2. All other set-up conditions are the same as in 
the case of the uniform mesh. Figures 6.8(a)-6.8(c) show the elevation plots at 
various times obtained by the adaptive implicit-explicit two-pass scheme (AIE-2). 
Figures 6.9(a)-6.9(c) show the distribution of the implicit elements for the AIE-2 
computations. The peak amplitude value at t = 0.72 is 0.969 for the implicit, EBE and 
AIE-2 schemes. The results for the explicit two pass (EXP-2) scheme are shown in 
Figures 6.10(a)-6.10(b). It can be seen in Figure 6.10(b) that in the left region of the 
domain the solution becomes unstable due to the high Courant number. 

The mean bandwidth (averaged over the number of time steps) of the global 
coefficient matrix for the AIE scheme in the case of the uniform mesh is 25% of that 
for the implicit scheme. The average number of implicit elements at each time step is 
173. For the nonuniform mesh the average mean bandwidth for the AIE scheme is 
79%, and the average number of implicit elements is 1036. 

Two-Dimensional Rigid Body Rotation of a Cosine Hill 
( Rotating Puffl 

The set-up conditions for this problem are the same as in the uniform mesh case 
of the translating puff problem except that the velocity field is rotational with respect to 
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Fig. 6.7 Elevation plots for the translating puff on a uniform mesh obtained by the 
EXP-1 scheme, (a) at t = 0.3. (b) at t = 0.72. 
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Fig. 6.8 Elevation plots for the translating 
puff on a nonuniform mesh obtained by 
the AIE-2 scheme, (a) Initial condition, 
(b) at t = 0.3. (c) at t = 0.72. 


Fig. 6.9 Distribution of the implicit elements 
for the AIE-2 calculations of the translating 
puff on a nonuniform mesh, (a) at t = 0.0. 
(b) at t= 0.3. (c) at t = 0.72. 
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Fig. 6.10 Elevation plots for the translating puff on a nonunifoim mesh obtained by the 
EXP-2 scheme, (a) at t = 0.3. (b) at t = 0.72. 
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the center of the domain (i.e., u t * - Xj + 0.5, u 2 = Xj - 0.5), and all boundary 
conditions are of the Dirichlet type and are homogeneous. The time step is adjusted to 
give a Courant number of 0.216 at the tip of the cosine hill. A full revolution is 
achieved in 200 time steps. Figures 6.1 l(a)-6.1 1(c) show the results at various times 
obtained by AIE-1. Figures 6.12(a)-6.12(c) show the distribution of the implicit 
elements at the corresponding times. The peak amplitude after a full revolution is 
found to be 0.984 for the implicit and EBE schemes, and 0.980 for AIE-1. The mean 
bandwidth (averaged over the number of time steps) for the AIE scheme is 26% of the 
implicit scheme. The average number of implicit elements at each time step is 173. 
As in the case of the translating puff (uniform mesh), the results for EXP-1 are not 
satisfactory. Figures 6.13(a)-6.13(b) demonstrate the performance of EXP-1 in this 
problem. 


Shock Structure/ Entropy Condition Test Problem 

The problem is governed by equation (2.1 1). A computational domain of unit 
length is taken. The number of elements is 40. The initial condition is shown in frame 
0 of Figure 6.14. The explicit solution is very close to the implicit solution such that 
there is no need to test the AIE scheme in this problem. EBE and the implicit schemes 
are compared. The differences are indistinguishable. The results obtained by the EBE 
scheme are shown in Figure 6. 14. Each frame corresponds to a time step of At = 0.6. 
As can be observed, the initial condition contains two stable shocks and one unstable 
shock. The unstable shock immediately breaks down, and eventually the two stable 
shocks merge to form a steady shock. Figure 6.15 shows the history of the EBE 
iterations (number of pseudo-time iterations in the first corrector step versus the 
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Fig. 6.1 1 Elevation plots for the rotating 
puff obtained by the AIE-1 scheme. 

(a) Initial condition, (b) at t = 3.14 . 

(c) at t = 6.28. 


Fig. 6.12 Distribution of the implicit 
elements for the AIE-1 calculations of the 
rotating puff, (a) at t = 0 . (b) at t = 3.14. 
(c) at t = 6.28. 
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Fig. 6.13 Elevation plots for the rotating puff obtained by the EXP-1 scheme, 
(a) at t = 3.14. (b) at t = 6.28. 
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Fig. 6.15 Shock structure /entropy condition test problem: History of the number of 
EBE iterations. 
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corresponding time). 


Flow Past a Circular Cylinder 

A two-dimensional viscous flow past a circular cylinder has been of special 
interests among researchers in the field of computational fluid dynamics ( see, e.g., 
[31,32] ). In this problem we have 1940 elements and 2037 nodal points. A refined 
and implicit zone is located around the cylinder. Diffusivity is 0.0025 giving a 
Reynolds number of 100 based on the diameter of the cylinder. The time step is fixed 
at 1.0. The number of corrections at each time step is allowed to be as many as the 
convergence criterion (10' 6 ) dictates.The mesh, the boundary conditions, and the 
distribution of the implicit elements are shown in Figure 6.16. Figures 
6.17(a)-6.17(c) show the streamlines, relative streamlines and isovorticity lines for the 
symmetric solution at time = 700 obtained by the implicit-explicit (IMEX) scheme. 
After that, an artificial disturbance is placed shortly to initiate a nonsymmetric solution. 
Figures 6.18(a)-6.18(c) show the nonsymmetric results at time = 1,200 obtained by 
the GBE scheme. The implicit and the implicit-explicit schemes all give very close 
results (differences less than 0.5%). The explicit scheme diverges. The mean 
bandwidth of the implicit-explicit scheme is 47% of the implicit scheme. The number 
of implicit elements is 320. 


Driven Cavity Flow 

Because of the discontinuity of the velocity at the comers, driven-cavity flow 
problem has become an important example for testing a numerical method during the 
recent years ( see, e.g., [33] ). In this problem, the flow velocity is zero on three 
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T = 1 ,co=0 


VP = 0.125 
© =0. 



am = o, 
a©/an = o. 


¥=-l ,©=0 


* ¥ = 0, ami =0 on the cylinder. 


Fig. 6.16 Flow past a circular cylinder at Reynolds number 100: Finite element mesh, 
boundary conditions and the distribution of the implicit elements for the IMEX 
calculations. 
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Fig. 6.17 Flow past a circular cylinder at Reynolds number 100: Symmetric solutions 
obtained by the EBE scheme at t = 700 (a) Local streamlines, (b) Relative streamlines, 
(c) Isovorticity lines. 
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Fig. 6.18. Flow past a circular cylinder at Reynolds number 100. Nonsymmetric 
solutions obtained by the EBE scheme at t = 1200. (a) Local streamlines. 

(b) Relative streamlines, (c) Isovorticity lines. 




sides and unity on the fourth side. A uniform mesh of 30 x 30 elements in a 1 x 1 
domain is chosen. Diffusivity is 0.002S, time step is 0.1, and Reynolds number is 400 
based on the side length of the square domain. The number of corrections at each time 
step is limited to 5. Figure 6.19(a) shows the mesh, the boundary conditions, and the 
distribution of the implicit elements. Results obtained by the IMEX scheme are shown 
in Figures 6.19(b)-6.19(d). Those results are indistinguishable from the implicit and 
EBE solutions (differences less than 0.001%) . The mean bandwidth for 
implicit-explicit scheme is 44% of the implicit scheme. The number of implicit 
elements is 403. 
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CHAPTER 7 


CONCLUSIONS 

In this report, we have presented approximate solution schemes for large 
equation systems resulting from finite-element formulation of fluid dynamics 
problems. The element-by-element (EBE) approximate factorization scheme is 
essentially an iterative scheme which totally eliminates the need for the formation, 
storage, and inversion of a large global matrix. All computations are performed at the 
element level. The method keeps the versatility of the finite-element formulation in its 
easy adjustment to irregular mesh and its easy implementation to boundary conditions 
and source terms. Futhermore, the EBE approximate factorization procedure is 
parallelizable and this aspect of it makes it especially favorable since parallel 
computations are expected to play a more and more important role in computational 
fluid dynamics. 

Implicit-explicit schemes in nature are approximations to implicit schemes, yet 
they substantially reduce the cost of formation, storage, and inversion of a large global 
matrix. In this approach, for the elements which are designated to be implicit, the 
element level matrices are kept as they are, whereas for the explicit elements the 
element level matrices are approximated by a diagonal matrix. In the adaptive 
implicit-explicit (AIE) scheme, the implicit elements are selected adaptively. Selection 
of the implicit and explicit elements is based on several criteria including the stability 
and accuracy. For the stability consideration, the elements for which the Courant 
number exceed the stability limit of the explicit scheme are assigned to be implicit. 
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while explicit elements are introduced elsewhere to reduce cost in computer memory. 
The accuracy criterion is determined by the dimensionless wave number of the solution 
from the previous time step ( or iteration). By appropriately numbering the nodes 
based on the distribution of the implicit elements, one can substantially reduce the 
bandwidth of global matrices. This scheme allows us to have implicit refinement 
where it is needed. Compared to other adaptive concepts such as element 
redistribution or element subdividing, the A1E scheme is far easier to implement 

We have applied these schemes to various problems governed by the 
convection-diffusion equation and the vorticity stream-function form of the 
two-dimensional Navier-Stokes equations. Streamline-Upwind/ Petrov-Galerkin 
formulations are employed for the spatial discretization of our model problems and a 
predictor-corrector algorithm is used to solve the resulting semi-discrete equation 
system. The results obtained by EBE, IMEX and AIE in all cases are 
indistinguishable from those obtained by the true implicit formulations while the 
storage required is substantially reduced. The savings in the storage of the coefficient 
matrix in various test problems are summarized in Table 7.1. Recently the author has 
developed a new finite-element procedure for vorticity stream-function formulation in 
which the coefficient matrix derived from the Poisson’s equation is symmetric and 
positive-definite [34]. For procedures which involve symmetric positive-definite 
matrices, convergence of the EBE, IMEX and AIE schemes are more predictable. 

The results indicate a great potential for the future work in this area. It is 
believed that it will not be very long before these schemes are accepted as powerful 
tools in large-scale computing. 
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Case 


Savings in the storage of 
the coefficient matrix 


Scheme 


Translating Puff 
(uniform mesh) 


Translating Puff 
(nonuniform mesh) 


Rotating Puff 


Flov past a Circular 
Cylinder 


Driven Cavity Flov 


75 % 


21 % 


74 % 


53% 


56% 


AIE-1 


AIE-2 


AIE-1 


IMEX 


IMEX 


Table 7. 1 Savings in the storage of the coefficient matrix for various test problems. 
* no renumbering of the nodes ( see Remark (2) in Chapter 5) 
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